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Abstract. - We present a mean-field dynamical theory for single semiflexible polymers which can 
precisely capture, without fitting parameters, recent fluorescence correlation spectroscopy results 
on single monomer kinetics of DNA strands in solution. Our approach works globally, covering 
three decades of strand length and five decades of time: it includes the complex cross-overs occur- 
ring between stiffness-dominated and flexible bending modes, along with larger-scale rotational 
and center-of-mass motion. The accuracy of the theory stems in part from long-range hydrody- 
namic coupling between the monomers, which makes a mean-field description more realistic. Its 
validity extends even to short, stiff fragments, where we also test the theory through Brownian 
hydrodynamics simulations. 



For a polymer in solution, hydrodynamics introduces 
long-range coupling between different points on the chain 
contour, its strength falling off with inverse distance like 
1/r. Though this has long been recognized as a crucial fac- 
' tor in understanding the dynamics of flexible polymers [1] , 
its importance in the case of semiflexible and stiff chains is 
not fully appreciated. The cause is often expediency: even 
without taking long-range coupling into account, the na- 
ture of the local interactions — governed by bending stiff- 
ness and inextensibility — presents a formidable problem 
in constructing a theory of semiflexible polymer dynam- 
ics. The most widely used model — the worm-like chain 
(WLC) [2] — yields nonlinear equations of motion, and thus 
all dynamical theories of the WLC have been approxi- 
mate. One approach is the weakly-bending assumption: 
deriving equations of motion from a perturbation analysis 
around the rigid rod limit [3-7]. This is particularly rel- 
evant for certain biopolymers like actin, where the large 
persistence length, l p ~ 0(1 /jm), means that a broad dy- 
namical regime, consisting of motion on length scales much 
smaller than l p , will be dominated by the bending stiff- 
ness. Yet for less rigid cases like double-stranded DNA, 
where l p m 50 nm, many empirical situations will involve 
complex cross-overs between stiffness-dominated and flex- 
ible regimes at different time scales. Weakly-bending ap- 
proaches cannot provide an accurate description of these 
cross-overs, which hydrodynamics makes even more chal- 
lenging to model: the fluctuation modes at all length scales 
are coupled due to the long-range interactions. 



The need for a comprehensive, quantitatively accurate 
theory of semiflexible polymer dynamics is made more 
urgent by advances in single-molecule experimental tech- 
niques. Fluorescence correlation spectroscopy (FCS) can 
already probe double-stranded DNA kinetics at the level 
of a single monomer [8-10]. The most recent results, by 
Petrov et. al. [10], reveal a rich sequence of dynamical be- 
haviors for the motion of a tagged end in DNA fragments 
varying in length from L m 30 — 7000 nm. In our work, we 
show that a single theory, without any fitting parameters, 
can give an excellent description of these experimental re- 
sults over the entire range, covering the whole cross-over 
between stiff and flexible chain dynamics. Our approach 
is based on a Gaussian mean- field theory (MFT) [11-15], 
which is in itself surprising: Gaussian models are usu- 
ally considered tools for flexible polymers, with limited 
applicability as one nears the rigid rod limit. What we 
demonstrate is that hydrodynamics — the complicating el- 
ement in the theory — is what underlies its success: the 
MFT becomes more accurate because of the long-range 
interactions. We specifically concentrate on highly stiff 
chains, the most difficult terrain for a Gaussian theory, and 
investigate the strengths and limitations of the MFT ap- 
proximation. Since the experimental results for short-time 
dynamics of stiff fragments are not well resolved, we ad- 
ditionally validate the theory through Brownian hydrody- 
namics (BD) simulations. Though the MFT is restricted 
to quantities which are spatially-averaged over all coordi- 
nate directions (i.e. the mean square displacements mea- 
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sured through FCS), it can capture the full complexity 
of the monomer motion: hydrodynamic effects, rotational 
and center-of-mass diffusion, and the cross-overs between 
dynamical regimes at different time scales. 

We begin by reviewing the general MFT approach to 
scmiflcxiblc polymer dynamics (a more detailed treat- 
ment can be found in Ref. [15]). Polymer stiffness is 
typically modeled by the WLC Hamiltonian, C/wlc = 
^IpksT j ds(d s u(s)) 2 , describing the bending energy of 
a polymer contour r(s), < s < L, with persistence 
length l p and tangent vector u(s) = <9 s r(s). The local 
inextcnsibility of the contour is expressed through the 
restriction u(s)| = 1 at each s. The partition func- 
tion Z is a path integral over all possible contours, Z = 
J Vu Y\ s <5(|u(s) — 1) exp(— /3CAvlc), with the S functions 
enforcing the inextensibility constraint and (5 = l/(fcsT). 
Since this constraint is nonlinear, only certain equilib- 
rium properties can be calculated exactly, and deriving 
any dynamical quantities requires an approximation. In 
our case, we obtain a more tractable Gaussian mean- 
field model by estimating Z through the stationary phase 
approach [12], Z ps exp(— (5Tmf) = / exp(-/?[/ MF ), 
where the MFT Hamiltonian Umf has the form: Umf = 
(e/2)Jds (&u(s)) 2 + V Jdsu 2 (s) + v (u 2 (Q) + u 2 (L)). 
Here local inextensibility has been relaxed, and the pa- 
rameters v and vo are related through the stationary 
phase co nditi on, 8 v Tmf = O^Fmf = 0. The latter 
yields \Jve/2 = vq = iksT/A. This is equivalent to 
making v and uq act as Lagrange multipliers enforcing 
the global and end-point constraints J ds(u 2 (s)) = L, 
(u 2 (0)) = (u 2 (L)) = 1. By setting the bending mod- 
ulus e = (Z/2)l p kBT, the Hamiltonian Umf can be 
tuned to reproduce exactly the tangent-tangent correla- 
tion (u(s) • u(s')) = exp(— \s — s'\/l p ) of a WLC with 
persistence length l p , as well as related thermodynamic 
averages like the mean square end-to-end vector (R 2 ). 

The Gaussian approximation provides a starting point 
for deriving the dynamics of the system, following a hy- 
drodynamic pre-avcraging approach similar to that of 
the Zimm model [13, 15]. The behavior of the chain 
contour r(s, t) obeys the Langevin equation, dtr(s,t) = 
- J ds' /iav g (s - s')SUuF/Sr(s',t) + £(s,t), where £(s,f) 
are Gaussian stochastic velocities whose components have 
correlations governed by the fluctuation-dissipation theo- 
rem: (&(M)&(s / ,t , )> =2kBT5 ij 6(t-t') f j lavs (s-s'). Here 
A*avg(s — s') is the pre-averaged mobility tensor, obtained 
from the continuum Rotne-Prager tensor '//(s, s'; x) [13, 
16] describing long-range hydrodynamic interactions be- 
tween two points s, s' on the contour at spatial separation 
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1 is the 3x3 identity matrix, a is a microscopic length 



scale corresponding to the monomer radius, 77 is the vis- 
cosity of water, (1q = 1/6^770 is the Stokes mobility of a 
sphere of radius a, and the Q step function excludes un- 
physical situations involving overlap between monomers. 
If a polymer configuration with points s and s' separated 
by x has an equilibrium probability G(s,s';x), then the 
pre-averaged mobility is defined through the integration: 
J d 3 x < jt(s,s'- 1 x.)G(s,s';x) — ^ a vg(s — s') 1 . In order to 
determine the relative importance of hydrodynamic ef- 
fects, we will also compare the free-draining case where 
the long-range interactions are turned off, and the mobil- 
ity reads /-t^ g ( s - s ') = 2a/J, S(s - s'). 

The Langevin equation can be solved through normal 
mode decomposition, yielding a set of coupled stochastic 
partial differential equations for the normal mode ampli- 
tudes. These can be diagonalized, resulting in a set of de- 
coupled normal modes characterized by relaxation times 
r„ (in decreasing order such that T\ is largest). As in 
Ref. [15], we impose a mode number cutoff M = L/8a to 
agree with BD simulations at short times, approximately 
modeling the discrete nature of the chain at distances on 
the order of the monomer radius a. Results at larger 
length scales (i.e. the experimental comparisons discussed 
below) do not depend on the details of the cutoff. 

To test the accuracy of the MFT, we compared the 
theoretical results to BD simulations [17] of a bead- 
spring worm-like chain consisting of N monomers of ra- 
dius a (contour length L = 2aN). The bead po- 
sitions Vi(t) obey a discrete Langevin equation, with 
monomers coupled hydrodynamically through the Rotne- 
Prager tensor [16]. The elastic potential of the chain 
U = Ubcn + U s ti + J/lj has three parts: (i) a bend- 
ing energy J7 bo n = {l p k B T/2a)J2i( 1 ~ cos9i), where B l 
is the angle between adjacent bonds; (ii) a harmonic 
stretching term C/ Str = (j/Aa) J2i ( r ;+M — 2a) , where the 
large modulus 7 = 200ksT/a ensures inextensibility and 
r j — Tj ; (iii) a truncated Lennard- Jones interaction 
w£i<,-e(2o - r id )[(2a/nj) 12 - 2(2a/r iJ ) 6 + 1] 
In our simulations we set N = 50, 
with various l p = 10a — 200a, and used a Langevin time 
step t — 3 x 10~ 4 a 2 / (kBT/io). Thermodynamic averages 
were derived from 5-20 independent runs, with data collec- 
tion occurring every 10 2 — 10 3 steps, and each run lasting 
10 8 - 10 9 steps. 

We start our analysis by focusing on one specific dynam- 
ical quantity: the mean square displacement (MSD) of a 
chain end-point, A ond (<) = ((r(L,t) — r(L, 0)) 2 ). The top 
panel of Fig. 1 shows A en( j(£) for DNA strands of various 
length, L w 100 — 20000 bp, taken from the FCS measure- 
ments of Petrov et. al. [10]. Superimposed is the MFT 
prediction, without any fitting parameters. The constants 
in the theory are taken from the experimental conditions 
and the literature: T = 298 K, 77 = 0.891 mPa • s, a = 1 
nm, a rise per bp of 0.34 nm, l p = 50 nm. The agreement 
between the MFT and experimental results is remarkable: 
the average discrepancy in the time range t = 10 _1 — 10 2 
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Fig. 1: Top: A en< i(t), the mean square displacement of an end- 
monomer in a dsDNA strand, for various lengths L = 98 — 
19941 bp. Bottom: the local slope a cn< j(t) = din A en d/dlnt of 
the log-log curves in the top panel. In both panels the circles 
are from the experimental FCS measurements in Ref. [10] . The 
solid lines are the MFT predictions, without any fitting param- 
eters: all values are taken from the literature and experimental 
conditions (see text). 

ms, where there is the least scatter in the FCS data, varies 
between 6 — 25% for the different L. This close agreement 
without fitting parameters is only possible if the full set 
of equations for the MFT normal mode amplitudes, in- 
cluding the off-diagonal coupling between normal modes 
due to hydrodynamics, is diagonalized and solved. If the 
off-diagonal elements are assumed negligible, as was done 
in Ref. [10], additional fitting parameters are required to 
get agreement: rescaling factors for the relaxation times 
and the diffusion constant. This in itself is a testament 
to the importance of hydrodynamic effects. In Ref. [15], 
we had compared the MFT to FCS data from a similar 
experiment [9] for a smaller set of DNA chain lengths: 
the discrepancies we observed between MFT and experi- 
ment at intermediate times are completely absent in the 
more recent and extensive FCS results analyzed here. This 
might indicate that the experimental setup or analysis in 
Ref. [9] may need to be re-examined. 

The local slope of the A en( j(t) curves in the log-log plot, 
(*end(£)j is shown in the bottom panel of Fig 1. a en d(t) is 
calculated for each time t by fitting a straight line to the 
log-log plot of data points within a small range of times 
ti defined by | log 10 ti/t\ < 0.15. The local slope would be 
constant for pure power law behavior, and what we find as 



L is increased is the gradual emergence of a scaling regime 
at intermediate times with a C nd ~ 0.62. The standard ex- 
pectation for long, flexible polymers in solution is given 
by the Zimm result, a on( j = 2/3 [1]. In our case, the sub- 
Zimm scaling and pronounced variation in a en d(t) with t 
is evidence that slow cross-over effects have a significant 
role [15]: the resulting deviation from the classical scaling 
theory of flexible polymers is thus experimentally observ- 
able, and precisely captured by the MFT. 

The time at which the slope reaches its minimum value 
is the same order of magnitude as n, the longest relax- 
ation time. For t > t\ we see a cross-over to center-of- 
mass motion, with a slope a C nd = 1- The oscillations in 
the FCS slopes in this cross-over are due to uncertainties 
in A enc j(t) at large t, arising from noise in the long-time 
exponential tails of the FCS correlation functions. There 
should be another cross-over at short times, where the 
length scale of the fluctuations is comparable to the per- 
sistence length, A ond (t) < l p « 2.5 x 10~ 3 /mi 2 . Here 
the expected behavior is a power-law scaling with expo- 
nent a en d = 3/4 in the free-draining case; hydrodynamics 
introduces logarithmic corrections to this scaling, which 
are observable as an increase in the exponent on the order 
of 10% [3,5]. Unfortunately in the time range where we 
should see this stiffness-dominated regime the FCS data is 
not well-resolved: there is too much scatter in the A en d(i) 
results when t < 10 _1 ms for us to be able to extract ac- 
curate local slopes. Note that the MFT curves are still in 
excellent agreement even for the shortest chains examined, 
where L is smaller or comparable to l p ss 50 nm (w 150 
bp). But what is being measured for these short, stiff frag- 
ments is essentially only the cross-over to center-of-mass 
diffusion, with the slope a on d approaching 1. 

We thus need an alternative approach to validate the 
MFT description of polymer dynamics in the cross-over 
to the stiff regime. In this case BD simulations are an 
ideal tool: computational constraints restrict us to rela- 
tively short chains, but this is precisely the limit we want 
to investigate, where l p ^> a, L ~ 0(l p ). Fig. 2 shows 
the comparison between simulation and MFT A cn( j(t) for 
a chain of length L = 100a and varying l p , going from the 
flexible case of l p /L = 0.1, to the stiff one of l p /L = 2.0. 
The regular theory /simulation results with hydrodynamic 
interactions (shown in blue) are contrasted to the free- 
draining case (shown in red). We see a quite interesting 
trend: the hydrodynamic MFT continues to be an accu- 
rate predictor of A en( j(<), even as l p becomes larger than L. 
In the time range t = 10 1 — 10 4 a 2 /ksT^Q the average error 
between the hydrodynamic MFT and simulation A cn d(i) 
varies from 3 — 15% for the different l p , similar to the er- 
rors seen in the experimental comparison above. On the 
other hand, the free-draining MFT is much less accurate: 
it noticeably overestimates the short-time A on( j(i), and the 
average errors compared to the BD results are 4-8 times 
larger than in the hydrodynamic counterparts. This is also 
plainly seen in the local slopes plotted in the bottom pan- 
els of Fig. 2: the non- hydrodynamic MFT performs signif- 
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Fig. 2: Top: A en d(i), the MSD of an end- monomer in a 
semiflexible polymer, for L = 100a and three different ra- 
tios lp/L (a = monomer radius). Bottom: the local slope 
Qend(^) = din A en d /dint of the log- log curves in the top pan- 
els. In all panels solid lines are MFT predictions, while circles 
are BD simulation results. The upper curves (blue) include 
long-range hydrodynamic interactions, in contrast to the lower 
curves (red) which correspond to the free-draining limit. 



icantly worse. Though hydrodynamics introduces another 
level of complexity into our approach, requiring an ad- 
ditional approximation in the form of pre-averaging, the 
resulting MFT is quantitatively more successful than the 
simpler free-draining theory This reflects a general well- 
known feature of mean- field theories: they are closer to 
reality in systems with long-range interactions. Due to hy- 
drodynamics, every point on the chain is coupled to every 
other point through the Rotnc-Prager tensor. For a free- 
draining polymer, the dynamics of a point on the chain 
are determined solely by the local bending/extensibility 
interactions with its nearest neighbors. Thus we can ex- 
pect a Gaussian mean-field description in this case to be 
a much cruder estimate. 

For the BD simulation results in Fig. 2, the evolution 
of the dynamics from the flexible to stiff limits roughly 
agrees with earlier scaling theory expectations: for the 
free-draining chain with l p /L = 0.1 we see an intermedi- 
ate time regime with a cn d approaching the Rouse value of 
1/2, appropriate for flexible polymers [1], which increases 
to around 3/4 as we move to larger l p . Adding hydro- 
dynamics has several effects: the center-of-mass diffusion 
becomes faster (due to hydrodynamic cntrainment) , the 
relaxation times become shorter, and the local slopes are 
all shifted upwards. The a on d plateau near 3/4 in the 
stiff limit is now near 0.85. For both the hydrodynamic 



and free-draining cases, note that only at the largest l p 
do we see the emergence of nearly pure power-law scal- 
ing behavior at short times: cross-over effects dominate in 
the more flexible cases. Moreover the shape of the a en d(t) 
curve changes as l p is increased: two local minima are 
formed, a shallow one in the short-time plateau region, 
and a much sharper dip at longer times on the order of T\ . 
This sharp dip, connected to the appearance of a new dy- 
namical regime, will be discussed further below. First we 
turn to a more basic question: given that at l p /L = 2.0 
we are approaching the stiff rod limit, how is it that a 
Gaussian model can still work so well? 

To answer this question, we look at the dynamics of 
semiflexible and stiff chains in more detail, by decom- 
posing the motion into components approximately par- 
allel and perpendicular to the chain contour. For the 
chain end-point, we define a tangent vector u avg (i) = 
(u(L,t) + u(£,0))/|u(£,i) + u(i,0)|, which is an av- 
erage between the end-point tangent at times t and 0. 
The parallel and perpendicular components of A en( j(i) are 
then given by: Aj) nd (t) = {[(r(L,t) - r(L,0)) • u a v g (t)] 2 ), 
A^j(i) = (A cnd (t) - A|) nd (t))/2. This decomposition is 
applied to BD simulation results of a chain with L = 100a 
and lp/L = 0.4, 2.0 in the top panels of Figs. 3-4. The to- 
tal MSD, A en d(i), and the corresponding MFT estimate is 
also plotted for comparison, and all quantities are shown 
both in the hydrodynamic and free-draining cases. At 
times greater than the largest relaxation time, t > t\, 
where the orientation of u(L,t) and u(L,0) are uncorre- 
cted, the components converge to the same value, equal 
to 1/3 the total A en d(i). For shorter times we see an 
anisotropy, becoming more prominent as the chain be- 
comes stiffer: AjJ nd (f) < A^ d (i), due to the suppression of 
parallel fluctuations with increasing l p . In the absence of 
hydrodynamics, one can derive scaling predictions for the 
two components based on the weakly-bending approxima- 
tion: A^ nd (t) should have an exponent a^ d = 3/4 [3,5], 
while the effects of tension propagation along the backbone 
lead to a higher exponent cr_ d = 7/8 for A| nd (i) [6,7]. 
We see these two dynamical behaviors most clearly in the 
free-draining BD results for l p /L = 2.0, with the local 

slopes a^ nd (t) and a|j nd (t) approaching 3/4 and 7/8 at in- 
termediate times (with slight modifications due to cross- 
over effects). Adding hydrodynamics changes the scaling, 
shifting both of these exponents up by 5 — 15%. 

The MFT entirely misses the anisotropy: given the 
isotropic nature of the Gaussian Hamiltonian, it predicts 
Ajj nd (i) = A^ d (i). However, the MFT will still work for 
any quantity which is averaged over all spatial directions, 
and thus can successfully estimate the total A cnd (£). In 
the stiff limit, the spatially averaged MSD is dominated 
by large fluctuations perpendicular to the backbone, for 
which the Gaussian MFT gives a reasonable description; 
in fact, as we will argue below, the MFT normal modes 
at large l p for n > 1 effectively behave like the trans- 
verse modes in the standard weakly-bending perturbation 
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Fig. 3: Top panels: the end-point MSD A cn d(t) and local slopes 
c*end(i) for a chain with L = 100a and l p /L = 0.4. Both hydro- 
dynamic and free-draining results are shown. MFT estimates 
for the total A en d(i) are drawn as solid lines. The BD simu- 
lation results show both the total A en d(i) and the decompo- 
sition into components, Aj nd (i), A^dM- Vertical lines mark 
the two largest relaxation times, n and T2, calculated from 
the MFT. Bottom panels: analogous results for the end-to-end 
MSD Aee(i) and corresponding local slope a cc (t). The pro- 
jected length MSD, A pro j = ((R(t) - R(0)) 2 } , is plotted as well 
for comparison. 



analysis. In contrast, the longitudinal fluctuations of a 
nearly stiff rod cannot be approximated well by a Gaus- 
sian model. But because their contribution to the to- 
tal A en( j(i) is small, the MFT manages to capture the 
spatially-averaged dynamics. While these limitations are 
not relevant to modeling the FCS results for A en d(i) of 
a freely diffusing chain, they would be significant for an 
anisotropic experimental setup, i.e. a chain under tension. 



Fig. 4: Same as in Fig. 3, but for l p /L = 2.0 



In this case the MFT would have to be refined [18]. 

Another way of looking at the chain internal kinet- 
ics is through the end-to-end vector MSD, A ee (t) = 
((R(i) - R(0)) 2 }, where R(i) = r(L,t) - r(0,t). For 
the same parameters used in the top panels of Figs. 3- 
4, the analogous A ee (i) results are shown in the bot- 
tom panels, together with the corresponding local slopes 
a ee (t). The decomposition here is defined with respect to 
the unit vector R avg (i) = (R(t) + R(0))/|R(i) + R(0)|, 
so that Al(t) = ([(R(t) - R(0)) • R a v g (i)] 2 ), A^(t) = 
(A ee (t) — Aee(t))/2. The component Al e (t) at short times 
is approximately equal to the MSD of the projected length, 
A proj (t) = ((R(t) - R(0)) 2 }, where R(t) = |R(f)|. This 
is a well-studied quantity for characterizing semiflexible 
polymer dynamics [3,19], and is also plotted in Figs. 3-4. 
In the stiff limit, the time scale at which A pro j(£) satu- 
rates is the upper bound for relaxation of both || and _L 
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wide range of flexibility, even for short, stiff fragments. 
The latter case is particularly interesting with respect to 
DNA: there are claims that the elastic energy of the WLC 
model may no longer be applicable at length scales < 100 
nm [21]. If the resolution of FCS experiments could be in- 
creased to probe this regime, it would provide a direct and 
independent test, together with our theory and simula- 
tion results, of DNA mechanical properties at these scales, 
highly relevant to cellular processes. 
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Fig. 5: Relaxation times T n versus mode number n calculated 
using the MFT for a semiflexible polymer with L = 100a and 
various l p /L. The dashed line shows an n~ 4 power-law scaling. 
Inset: same data, but plotted as r n l p versus n. 
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internal contour fluctuations. As can be seen in Figs. 3- 
4, this saturation time coincides approximately with t<i. 
Though A^(t) continues to grow between T2 and n, the 
dominant contribution in this range is rotational diffusion 
of the polymer backbone. The underlying reason is that 
for lp/L — > oo, the n = 1 normal mode in the MFT be- 
comes a purely rotational mode [14,20]. We can see this 
directly in the plot of r„ versus n for L = 100a and various 
lp/L in Fig. 5: as we move to stiffcr chains, n approaches 
a constant value, independent of l p . Up to a prefactor, 
this constant agrees with r r ~ L 3 / ^ksT ln(L/2a), the 
rotational relaxation time of a stiff rod of length L and di- 
ameter 2a. The behavior of r„ for n > 1 is quite different: 
the times scale approximately like ~ L i /l p n i (the inset of 
Fig. 5 shows data collapse for different l p ), with correc- 
tions due to hydrodynamics. This is exactly the predicted 
behavior of the relaxation times for transverse fluctuation 
modes in the weakly bending approach [3]. 

Thus we can now understand fully the various dynami- 
cal regimes seen in both the simulation and MFT results 
for A en d(i) in the stiff limit. For short times t < T2 the 
motion is dominated by the _L bending modes of the chain, 
and we see the characteristic 3/4 scaling (plus hydrody- 
namic corrections) in the total MSD curves. Between T2 
and t\ we have a regime controlled by rotational diffu- 
sion: the local slopes increase after T2 , dip sharply through 
T\, before rising to 1 at times t 3> t\, where center-of- 
mass translational diffusion is dominant. This rotational 
regime has not been correctly reproduced by any earlier 
theory based on the weakly-bending assumption. Though 
the MFT is limited to spatially-averaged properties, it 
does describe quantitatively the full cross-over between 
all three regimes, particularly important when making de- 
tailed comparisons to experiments. 

In summary, we have presented a mean-field approach 
to semiflexible polymer dynamics that gives a highly ac- 
curate description, without fitting parameters, of both ex- 
perimental FCS measurements on DNA and BD simula- 
tions. It incorporates hydrodynamics and works over a 
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